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A variational approach for the free energy is used to study the three-dimensional anisotropic XY 
model in the presence of a crystal field. The magnetization and the phase diagrams as a function of 
the parameters of the Hamiltonian are obtained. Some limiting results for isotropic XY and planar 
rotator models in two and three dimensions are analyzed and compared to previous results obtained 
from analytical approximations as well as from those obtained from more reliable approaches such 
as series expansion and Monte Carlo simulations. It is also shown that from this general variational 
approach some simple assumptions can drastically simplify the self-consistent implicit equations. 
The validity of the low temperature region of this approach is analyzed and compared to Monte 
Carlo results as well. 

PACS numbers: 75.10H 
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1. INTRODUCTION 

Ising and Heisenberg models are spin Hamiltonians 
which have been originally proposed in the study of mag- 
netism in pure and diluted materials [H, [3, d, Q ■ They 
are widely studied and employed in both classical and 
quantum contexts. However, a great interest has also 
been devoted to the XY model due to the richness of its 
phase diagram and the character of its transition Q . 
In particular, the XY model in two dimensions is of prime 
interest in the field of statistical mechanics because of its 
intriguing and unusual phase transition. Recently, the 
study of the XY model has been motivated not only 
due to its applicability in describing real magnetic ma- 
terials (for instance, it has been quite recently shown 
that the phase diagram of the compound RbFe(Mo04)2 
is in good agreement with the theoretical predictions for 
a two-dimensional XY classical model 0) but as well as 
in condensed matter systems, for example liquid crystals 
[§| and superconductors TO, 11|, among others. 

The one-dimensional version of the quantum spin-1/2 
XY model has been exactly solved by Lieb, Shultz and 
Mattis Although the quantum nature of real mag- 
nets cannot be forgotten, the study of classical models 
continues to be an important subject of research. The 
two-dimensional classical version of the XY model has 
been treated by Berezinskii and Kosterlitz and Thouless 
[1, |6|] and they showed that the vortices and anti- vortices 
play a central role in the thermodynamic behavior of the 
model. For higher dimensions this model (as well as a 
planar rotator version) has been studied through approx- 
imate analytical techniques 
Carlo simulations 



T = 
one 



is transformed into a classical d 



1 dimensional 



[Ij, IJJ, IJJ, IJJI and Monte 
17j |. In addition, classical models in 
general, and in particular the XY model, have become 
very popular recently in the context of quantum phase 
transitions, where a d-dimensional quantum model at 



The effort for a better theoretical understanding of 
the high temperature superconductors has also lead to 
an increase in treating anisotropic XY type models 
11, [23, 21, specially its two-dimensional version. 



However, despite the layered structure of these systems, 
the high temperature superconductors are not strictly 
two dimensional. For this reason, inter layer interaction 
should be important in describing their thermodynamic 
properties ll| , as well as the inclusion of possible crystal 
field interactions. 

In this work we study the classical anisotropic XY 
model with a crystal field interaction described by the 
following Hamiltonian 

H = -jJ2{S^~Sp+Spl} 

{f,r ) 

- J.Y. {^^^'^Z + SlSl} + DY,{S'f? , (1) 

where J is the exchange interaction between spins in the 
layers parallel to the xy plane and Jz the exchange in- 
teraction between spins in different adjacent layers. D is 
the crystal field and S" are the a = x,y,z components 
of a classical spin jS'pl = 1. The first sum runs on nearest 
neighbor spins {r,r) within the layers and the second 
sum on nearest neighbor spins (r , r ) between layers (in 
order to avoid confusion here and in the following ex- 
pressions, the subscript of the exchange interaction or 
the variational parameter in front of the sums will define 
whether the spins belong to the same layer or to adjacent 
layers). The last sum is made over the entire N spins on 
the simple cubic lattice. 

The Hamiltonian ^ can be written in a more con- 
venient form by means of a polar representation for the 
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spins 



Sp = (52 , S"! , S'5) = (sin 6'f^cos (l)p, sin 6*? sin (f>p, Sp), 



Sp = (^l-(5£)2cos#,^l-(55)2sin 0,-55) , (2) 

where 9p and 4>p are the spherical angles of the spin at 
the site f. In this representation the Hamiltonian given 
by Eq. ([1]) takes the form 

J 



^ = - 9 E V 1 - (S'^rJl - (55+,-)2 COs(</>.-+s 



Jz 



^E(^^-)'' 



(3) 



where a labels the four nearest-neighbor sites of f in the 
xy plane and c the two nearest-neighbor sites of r along 
the z direction. 

The above model has been treated according to the 
self-consistent harmonic approximation (SCHA) in the 
case D = for the three-dimensional model [l6| as well 
as the planar rotator model 15[ (in this case there is 
no z spin components). The approach in [l6| . however, 
only considers the cases Jz ~ J and Jz « J, that is, 
quasi-isotropic case or weak layer coupling and no general 
treatment has been done for any value of < Jz < J. 

In what concerns the theoretical approach, it is well 
known that the SCHA has been extensively used in the 
literature [ul, [3, 0, IISS H S S 111, 111, 0, H llil . 
It consists of replacing the original Hamiltonian ([3]) by a 
harmonic one by expanding the corresponding cosines to 
the second order with effective exchange constants that 
take into account the nonlinearities of the interactions. 
The effective couplings are then chosen to minimize the 
corresponding free energy of the system. In the present 
work we use a variational method based on Bogoliubov 
inequality for the free energy to study the model Hamil- 
tonian ID) in order to obtain its low temperature thermo- 
dynamic properties. This procedure is the same as the 
SCHA employed in previous works. Some simplifications 
are also suggested in order to easier handle the awkward 
self-consistent implicit equations. 

The plan of the paper is as follows. In the next sec- 
tion, we present the analytical procedure according to 
the variational approach for the free energy in order to 
obtain the thermodynamic properties of the model. In 
section 3, the numerical results are presented, where the 
role played by the anisotropy in the presence of a crystal 
field is shown to be for itself relevant (for instance, in 
the context of high temperature superconductors). Ad- 
ditional assumptions are also suggested which can make 
the implicit self-consistent equations easier to be han- 
dled. Some concluding remarks are discussed in section 
4. 



2. VARIATIONAL APPROACH FOR THE FREE 
ENERGY 

The present variational approach is based on the Bo- 
goliubov inequality for the free energy 



F<Fo + {H-Hoh))o = Hl), 



(4) 



where Ti. is the Hamiltonian in study given by Eq. ([T]) 
or 7^0(7) is a trial Hamiltonian which can be exactly 
solved and depends on variational parameters 7. F is 
the free energy of the system described by Ti, Fo is the 
free energy of the trial Hamiltonian Tio, and the thermal 
average < ... >o is taken over the ensemble defined by Ho- 
The approximate free energy is given by the minimum of 
$(7) with relation to 7, that is, F = ^min{j)- 

In general, the trial Hamiltonian should resemble, in 
some aspects, the one under study. In this case, Tio can 
be chosen as a sum of two parts 



Ho = Hi + HI, 



(5) 



in such a way that the first part is a kind of a planar 
Hamiltonian 



^^-iE(^^-+'^-'^^^)' + tE( 



and the second term is an axial Hamiltonian 



Hl = {D + 2J + Jz)Y,{Sif 



(6) 



(7) 



where 7 and 7^ stand for the variational parameters. 
This harmonic choice for Ho is also motivated by the fact 
that at low temperatures the angle differences \(j)r+a ^ 
4>p\ « 1 and \(j)p^g— « 1 so the cosines in Eq. ^ 
can be expanded up to second order to give the terms 
{4>f+a ~ fprY {4>r+c~ "^0- Sincc in this case 

there is a negligible variation in the z components of 
the spins, we can further assume that 55 « 55_|_^ and 
55 w 55_j^-, so the square roots in Eq. ^ are eliminated 
and we end up to a term D + 2J -\- Jz in the corresponding 
axial term. Thus, it turns out that this approach will be 
valid only in the low temperature region. 

Obtaining the averages regarding the trial Hamiltonian 
requires diagonalizing The planar Hamiltonian Hq 
can be diagonalized in the reciprocal space through the 
Fourier transform of (hp 



with the inverse transform (p^ given by 



(8) 



(9) 
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Using the transformations ([8]) in Eq. ([6]) and the fact 
that the system is translationally invariant we get 



q a 



(10) 



Summing now over the vectors a in the xy plane and c 
afong the z direction and rearranging terms we obtain 
the diagonal form of the planar trial Hamiltonian 



(11) 



where 7, = 7(2-cosgj;a-cosgya), 7^^ = 7^(l-cosg2c), 
a = |a|, c |c| and \(f>q\'^ = 4>q(t>-q- 

The diagonalizing of the axial term of the harmonic 
Hamiltonian ([7]) is obtained by introducing the corre- 
sponding Fourier transform of the z component of the 
spins Sp 



q 



— iq-r QZ 



and the inverse transform given by 



^q- ^2^^ ^r- 



Applying (HH) to Eq. ^ one gets 

ni = {D + 2j + j,)Y,\si\ 



(12) 



(13) 



(14) 



The mean value < T^o >o can be evaluated by using 
the equipartition theorem resulting in 

(^o)o = — ^ \ ^ — = NksT. (19) 

On the other hand, the mean value of < >o is not 
so straightforward computed. It can be written as 

J ■ 



E( V 1 - (^.-) - (^.5+a)'>o(cOS(0.-+S - 0f.))o 



J. 



r.c 



(20) 



where in the first two sums the mean value of the product 
of the planar and axial terms regarding the trial Hamil- 
tonian has been factorized and 



2{D + 2J+J,) 



(21) 



is the out-of-plane spin fluctuation. 

For the Gaussian variables (^f+s ~ (f'r) we can write 

(cos(0,-+s - #))o = e-^«"*^^+»— (22) 

where the mean value appearing in the exponential is 
given by 

((#+,- - cj^,r)o = I E(l - \-)(l^-|')o- (23) 



where {Sil"^ = 555^-. 

The partition function Zq can be computed by first 
noting that the planar and the axial parts of the harmonic 
Hamiltonian are independent so 

Zo = Tre-"^^' = Tre-^^^'+^i^ = Z^Z^. (15) 

Since both "Hq and Tig are quadratic in their variables 
one has 



where A, = ^{cos q^a + cos qya) and 



(24) 



q\ 1^ 0/ I ^ ■ 
2(79 +7?^) 

Similarly, for the Gaussian variable {(t>p+^ — 4>p) we have 



(cos((/.,-+5- (/.,-))o = e-5((*.^+c-0.)=>o^ 



and 



Zt = Tre 



n 



(i{lq+lqz) 



(16) 



{{4>r+c~4>rfh = ]^E(1 - ^|)(l'^d')0, (25) 



and 



(17) 



where Xq^ — cosg^c. In this way, Eq. (|20p assumes the 
form 



)(i-#',r>o 



where VI = D + 2J + Jz. The free energy Fq is then given 

by 



Fo 



kuT 



El- 



P{lq+lqz) 



J: 



r,c 



)(I0,|'>O 



(26) 
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where we have used an additional assumption that 5*5 « 
Spj^^ and Sp « Spj^g. The last term can again be com- 
puted from the equipartition theorem and as the terms 
in the sums do not depend on the respective indexes we 
have 

~2JN{1 — ((5'£)2^^^g-2Tv 5Z<r(2-(cosgj;a+cos(;„a))(|0,|^)o 



2v„ NDknT, 



-J,N{1 - ((55)2)o)e-^ i:,-(i-cosq.c)(|0,|^)o _^ 
The right hand side of Eq. ^ is then written as 



2n 



{27) 



■Vln^^ 



Vln — 



— 2JN{1 — ((55)^)o)e~ 27v E5(2-(cosg^a+cosg„a))(|0<,p)o 

Minimizing the above equation with respect to the 
variational parameters gives an upper bound limit for the 
free energy. The variational parameters are determined 
from the conditions 

^*^^'^^)^0(a) and = 0. (6) (29) 



97 



The mathematical expressions for Eqs. (|29|) are rather 
lengthy to be reproduced here. However, factorizing 
terms that can be canceled out and defining 

'yxy = {cos qxa + COS qyo)] (l^gHo, (30) 



(35) 



which are now two parametric equations from which the 
two variational parameters can be obtained. These equa- 
tions can be put in a more convenient form by noting that 



exp 



l-^E— 

I iV'Y2[27(l 



kBT[l - \) 



<36) 



Taking the continuum limit of the above 
equation in cylindrical coordinates — > 



2 c27r 

a c 



27r4 



limit (f « we find 



d9 I qdq / dqz and the long wave length 



exp ■ 



ksT 

" 27 



arctanga 1 , 1 . 

352 6 6 g 



(37) 



where g = ^ . Analogously we find 



exp{- 



■ 27 



3- 2arctan(72 -(_ (^2 _ g2 ln(l -| — 

5 3g2 V ■ 9 



},(38) 



so that the variational parameters are finally obtained 
from 



7 = J(l-((55)^)„)exp{ 



fcsT .arctan (72 
3g5 



(39) 



we arrive at the following expression 



[l-((5£)2)o] (2J,7. 



(31) 



(32) 



It is interesting to notice that this very same equation 
is obtained either from condition (|29b .) or ((29b ). mean- 
ing that both variational parameters cannot be obtained 
from this equation alone. However, from the Gaussian 
variable definitions (|22 p - ((25|) one can deduce the follow- 
ing additional relation for the fluctuations rjxy and rjz 



"^IVxy + IzVz 



kuT 



(33) 



Comparing now Eqs. ((32)) and ((33|) we end at the follow- 
ing identifications 



7 = J l-((55)^)o e 



(34) 



Iz = J.(l-((^5)")o)exp{ 



" 27 ^5 



1 



(2arctan55 -f 52 - g2 ln(l + -))]}. (40) 



The two equations above can also be obtained by em- 
ploying the usual SCHA. Thus, for a given value of 
t = keT/J, D/J and Jz/J one can solve the non-linear 
system ((39)) and (00)) to get 7/ J and 7z/ J, and from them 
the desired thermodynamics of the model. For example, 
by taking the continuum limit in the long wave length 
regime the x component of the magnetization is given by 



,^ , , fc^T ,arctan (02 ) 

™ = (1- o(('5F)')o)exp{- ^ ' ' 



2Tr^j' 



+ hnil + -)]}■ 
2 9 



(41) 



The transition temperature is obtained when the only 
solution of Eqs. ((23) and ((10)) is 7 = 7^ = 0. 
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TABLE I; Reduced transition temperatures tc for the 
anisotropic XY model in some isotropic limiting cases ac- 
cording to Monte Carlo (MC), series, present approach, and 
Assumption 1 and Assumption 2 (Ass, 1-2). 



KJK of Ref. ^ and 





MC or series 


present 


Ass. 1-2 


L> ^ oo, = 


0.90[32] 


1.472 


1.472 


D ^ oo, Jz — J 


2.17[33J 


2.190 


2.207 


D = 0, = 


0.78(2)[34] 


1.076 


1.076 


D = 0, Jz = J 


1.54(l)[16l, 1.55[35] 


1.605 


1.613 



Before exploring possible simpler solutions of Eqs. (IHS)) 
and (|40p and presenting the numerical results, it is worth- 
while now to discuss and compare some limiting cases. 

Let us consider first D ~ 0. In this case we have the 
anisotropic XY model. The out-of-plane spin fluctuation 
takes the form 



AJ + 2J^ 



(42) 



which is the same as that obtained by Costa et al. 
[l6| . The transition temperature for the two-dimensional 
model Jz = 0, as well as the isotropic three-dimensional 
model Jz = J, are also identical to those from reference 
[l^ . respectively, tc — 1.076 and = 1.605. However, 
it is worth to stress here that the present approach not 
only is a generalization over the model treated by Costa 
et al. [16l | to other values of the axial interaction J^, 
but also includes, as it will be seen below, the important 
contribution of the crystal field interaction D. 

In the limit D — s- oo, as it will be discussed be- 
low, we have the planar rotator model. Again, we find 
the same transition temperature for the two-dimensional 
model tc = 1.472 and, for the three-dimensional model, 
tc = 2.190 [l^. The transition temperatures for _D = 
and D ^ oo are given in Table U for the isotropic two- 
and three-dimensional models together with those com- 
ing from other approaches. 

In addition, the present approximation can also be 
used to evaluate the behavior of the x component of the 
magnetization for small axial interaction Jz J and any 
value of D. In this limit, we also have g ^ so Eq. (PT|) 
gives 



7 



(43) 



For the planar rotator model we have {{SZ)'^)o — and 
the above equation should be compared to 



Sir if 



(44) 



(45) 



obtained from spin wave theory where now 72/7 plays 
the role of Jz/J of ReL The exponent here is not 

the same, although numerically comparable to other ap- 
proaches. 



3. Numerical Results, Parametric Procedure and 
Simple Assumptions 

Numerical Results 

In Figure[T]we have the reduced transition temperature 
as a function of the anisotropy ry = Jz/J for D = 0. From 
hereon the present approach means the results obtained 
by solving the non-linear system of equations (j39p and 
(|40)) numerically (by using standard iterative procedures) 
for the variational parameters. As can be seen from this 
figure the results are quite similar to those from reference 
[l6| in the region close to the isotropic three-dimensional 
case Jz ^ J. On the other hand, the slope of the phase 
boundary close to isotropic two-dimensional case Jz = 
is zero according to the present approach while the 
procedure from reference [16| for the XY model (and 
also for the planar rotator model [l5| ) furnishes a positive 
slope. 




obtained from a SCHA where 7z/7 plays the role of 



0.4 0.6 
ri=J/J 

FIG. 1: Reduced transition temperature t — ksT/J as a 
function of ?7 = Jz/J for D = 0. The solid line represents the 
present approach and the dot-dashed line the results from Ref. 

Figure [2] shows the phase diagram in the reduced tem- 
perature t = ksT/J and rj = Jz/J plane for several 
values of D/J. One can clearly note that the crystalline 
anisotropy D plays an important role in the critical be- 
havior of the model. 
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For positive values of the crystal field D > the tran- 
sition temperature increases as D increases, since in this 
case the out-of-plane fluctuations are reduced implying a 
greater tendency of the spin components to lie in the xy 
plane. For D —^ oo we recover the planar rotator model 
because there will be no z component of any spin. As 
shown in Table HI in this limit, the result tc = 1.472 
for Jz — are quite different from that obtained by 
Monte Carlo simulations tc = 0.90, reflecting the fact 
that the present variational approach does not take into 
account vortices effects, which are relevant for such two- 
dimensional model. On the other hand, for — J the re- 
sult tc = 2.190 for the isotropic three-dimensional planar 
rotator model are quite comparable to the Monte Carlo 
simulations tc — 2.17. Moreover, in the two-dimensional 
limit (J^ = 0), Eqs. and (gO]) can be written as 



ksTc 



tc = 



4f 



2 + e{2+^y 



(46) 



which gives the transition temperature of Fig. [2]for rj = Q 
down to the value -j ~ —2. 




at which the out-of-plane fluctuation ||2T|) diverges. 

Figures [3] and 2] show the temperature dependence of 
the magnetization obtained from Eq. (|41[) where one 
can see a discontinuous behavior. This discontinuity in 
the magnetization is an artifact of the present method. 
The transition temperature is given when the non-linear 
system of equations ([M)) and (|^ admit only the triv- 
ial solution. In all cases they do not go smoothly to 
the trivial solution presenting thus a discontinuity. This 
is a general feature of such methods. For instance, the 
temperature behavior for D oo and Jz/J = 0.1 and 
Jz/ J — I shown in Figure 31 is quite similar to those given 
in Figure 3 of reference [l5|. 
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FIG. 3: Magnetization m as a function of the reduced tem- 
perature ksT/J for 77 — Jz/J = 0.1 and two different values 
of the crystal field, D/J = and D/J — 1, according to the 
present approach (full lines) , and Assumption 1 (dashed lines) 
and Assumption 2 (dot-dashed lines). 



Parametric Procedure and Simple Assumptions 



FIG. 2: Reduced transition temperature t — ksT/ J as a func- 
tion of rj = Jz/J for several values of D/J (indicated by the 
numbers) according to the present approach (full lines), and 
Assumption 1 (dashed lines) and Assumption 2 (dot-dashed 
lines) . 

For negative values of the crystal field D < we have 
an inverse situation. We note that as D decreases the 
transition temperature also decreases. This can be un- 
derstood because the out-of-plane fluctuations are now 
enhanced implying in a tendency for the spins to lie out 
of the xy plane and to become more Ising like. For a given 
value of D there is a critical value of the ratio rjc = Jz/J 
at which the temperature goes to zero. This value is 
given by 



Vc 



-D/J -2 



(47) 



If we take the ratio between the variational stiffness 
g = —, we can compute the quadratic fluctuations given 
by Eqs. pop and pip by taking the continuum limit in 
the long wavelength regime 



where 



Fig) 



27 



27 



1 - 2F(.g) 



arctan(g2) 1 1 l 



52 z z g 

On the other hand, Eq. ([5^ can be written as 

-t 



G(g) = Q In Q [2 -f 77aQ"-i] =^ 



(48) 



(49) 



(50) 
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1.5 



2.5 



FIG. 4: Magnetization m as a function of the reduced temper- 
ature ksT / J for D oo and two different values of 77 = Jz/J, 
7] = 0.1 and ri — 1, according to the present approach (full 
lines), and Assumption 1 (dashed lines) and Assumption 2 
(dot-dashed lines). 



where Q = e 



and ?7 = ^. From Eqs. ([48] 



-v^y ^ and ?7 = 
the parameter a and, in the same way from Eqs. 
and (j35p the quantity Q, can be written in terms of the 
ratio g as 



1 - 2F{g) 



(51) 



Finally, the temperature can be expressed in terms of 
G{Q) as 



t 



1 



-2G{Q) 



(52) 



So, knowing a priori the Hamiltonian parameters rj and 
D/J and for a given value of g one gets 



^(.9) 



G(Q) 



t. 



(53) 



This is indeed what really happens when we use the nu- 
merical solution of the implicit self-consistent equations 
([55)1 and (HOI). However, it would be quite nice if one 
could track the inverse path of Eq. ([53]) . Unfortunately, 
as one can see, this is in fact not possible since for a 
given t (the usual parametric procedure to get the dif- 
ferent thermal dependences) we cannot solve Eq. ([50]) 
because a is unknown. It is at this point where some ad- 
ditional assumptions could be made in order to explore 
simpler solutions of Eq. (|50|) . 



Assumption 1. a — 11 

In this case Eq. ((50|l reduces to 
QlnQ [2 + r]^Q"-^] 



2 [1 - ((^5)^)0] 



(54) 



In Figure [2] it is shown, by the dashed lines, the corre- 
sponding critical temperatures according to this assump- 
tion. It can be seen that, within that scale. Assump- 
tion 1 and the present variational approach are almost 
indistinguishable. For D = the transition temper- 
ature obtained from Eq. ([M]) for the two-dimensional 
model Jz = is identical to that from the present one 
tc = 1.076. However, for the isotropic three-dimensional 
model Jz ~ J one has tc — 1.613, which is slightly differ- 
ent from the value tc — 1.605 according to reference [3]. 
On the other hand, in the limit 13 —> 00 we find the same 
transition temperature for the two-dimensional planar ro- 
tator model tc — 1.472 and, for the three-dimensional 
case, tc = 2.207 is comparable to tc — 2.190 obtained 
from Rcf. 16]. These values are depicted in Table [J to- 
gether with those coming from other approaches. The 
magnetization as a function of the temperature is pre- 
sented in Figures [3] and [H A different behavior is ob- 
served, except for the planar rotator model. 

The quite good agreement of the present approach and 
Assumption 1 reflects the fact that the term -qaQ"^^^ in 
Eq. (|50p varies very weakly with g and thus with tem- 
perature. Reproducing the correct limits at 77 = and 
1 the phase boundaries are then a sort of smooth inter- 
polation for different values of the anisotropy, explaining 
the accidental good agreement. With this in mind other 
simplifications can also be done. 



Assumption 2. F{g) 



1 

2+g 



In the interval [0,1], the function F{g) given by Eq. 
(I49p can be approximated by (with an error of order less 
than 10%) 

1 



F{9) = 



9 



(55) 



This implies that a = 1 and rjz = rjxy, meaning that the 
quadratic fluctuations are not very different in vertical or 
in-plane bonds. The self-consistent Eq. ([50]) adopts now 
a very simple form 

Q'-QP+"'- 2[i-«U.] ' 

Unlike Eq. (j54p the above equation permits obtaining an 
analytical expression for the critical temperature. Since 



dQ 
dt 



at tc one gets Qc 



tc = 



— e and 

2(2 + ry) 

i{2 + v) + e 



(57) 



The above equation also yields exactly the same limiting 
results as in the last column of Table [J The boundaries 
as a function of D are given in Figure [2] by the dot- 
dashed lines. In this case, for positive values of the crystal 
anisotropy the agreement is not as good as for negative 
values. The corresponding magnetizations as a function 
of the temperature are also shown in Figures [3] and U) 
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4. CONCLUDING REMARKS 



The anisotropic XY model in a crystalline field has 
been studied according to a variational approach for the 
free energy. This system is a generalization of the planar 
rotator model and the aniso trop ic XY model previously 
treated in the literature [H, • We believe we have ob- 
tained a satisfactory picture of the thermodynamic be- 
havior of the model as a function of its parameters. Not 
only is the phase diagram in the anisotropic case Jz ^ J 
different from the previous one obtained for D = 0, but 
the crystal field has been show to play an important role 
in the critical behavior of the system. 

Regarding Assumptions 1 and 2, one could see that the 
former one (a = rf) would imply that the mean quadratic 
fluctuation on a given bond scales with the coupling con- 
stant of the bond, while the latter one (a = 1, -q = 1) 
means that the quadratic fluctuations are not very dif- 
ferent in vertical or in-plane bonds. So, according to 
Assumption 2, the average contribution to the internal 
energy from each bond scales approximately with the 
coupling constant of the considered bond. This seems 
more acceptable than the conceptually wrong assump- 
tion a = rj, since one expects that softer bonds should 
develop stronger fluctuations. Thus, at flrst sight, from 
the numerical point of view, the better agreement of As- 
sumption 1 with the complete approach can be ascribed 
to an artifact of the method. 

It is clear from the approximation employed in this 
work that it should be valid only at low temperatures. It 
is also rather surprising that the values for the transition 
temperatures shown in Table [J for the three-dimensional 
model are quite comparable to those coming from more 
reliable methods (for the two dimensional model one 
would not expect such agreement due to the vortices 
effects). It would then be quite nice to have a clearer 
picture of the range of the temperature validity of the 
present procedure. In Figure [5] we show the out-of-plane 
spin fluctuations given by Eq. (PT|) for several values of 
the Hamiltonian parameters compared to Monte Carlo 
simulations. We have chosen this quantity because it is 
easily obtained from the present approximation, Eq. pi|) . 
The simulations have been done in the three-dimensional 
model described by Hamiltonian ([1]) and employing the 
single spin-flip Metropolis algorithm. A flnite lattice of 
Lx Lx L {L = 14) has been used with periodic boundary 
conditions. The averages have been taken by considering 
10^ Monte Carlo steps (MCS) per spin after equilibra- 
tion. For the flrst temperature, lOOL^ MCS have been 
discarded and, with the previous conflguration being the 
initial conflguration for the next temperature, additional 
3 X 10'^ have been discarded. It is evident from this fig- 
ure that the transition temperatures are not far from the 
temperature range where the corresponding results are 
comparable to the Monte Carlo simulations. 
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FIG. 5: Out-of-plane spin fluctuation as a function of the 
temperature for several values of the Hamiltonian parame- 
ters. The full lines are the results from Eq. (I21II and the 
different symbols are Monte Carlo simulations [33] • The ar- 
rows indicate the transition temperatures obtained from the 
present approximation. The MC results saturates at 1/3 as 
T ^ oo. 
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